\(~\)

This notebook helps visualize and explore surface and groundwater observations stored in the National Water Information System (NWIS). It was written to perform automated pulls and interactive mapping and plotting of data as it was being collected and analyzed for a specific project, but it can be modified for any region where there is data in NWIS.

\(~\)

1. Map sites

# watershed screening sites 
screensites <- whatNWISsites(bBox = c(-70.448, 41.626, -70.392, 41.677))
nums <- "504|505|506|507"
screensites <- screensites %>%
  filter(grepl(nums, station_nm))

# local demonstration sites
sites <- whatNWISsites(bBox = c(-70.405, 41.669, -70.395, 41.677))

# site info
baseurl1 <- "https://nwis.waterdata.usgs.gov/usa/nwis/qwdata/?site_no="
baseurl2 <-  "&agency_cd=USGS"

siteinfo <- readNWISsite(sites$site_no) 
#comment(siteinfo)
siteinfo <- siteinfo %>%
  dplyr::select(site_no, station_nm, site_tp_cd, dec_lat_va, dec_long_va, dec_coord_datum_cd, alt_va, alt_datum_cd, well_depth_va)%>%
  mutate(id = substr(station_nm, 9, 11)) %>%
  mutate(weblink = paste(baseurl1,site_no,baseurl2,sep = " "))

# reclassify site types for mapping
siteinfo$site_tp_cd[grep("M01", siteinfo$station_nm)] <- "multilevel sampler (MLS)"
clust <- "505-0|517-0|524-0|525-0|526-0"
siteinfo$site_tp_cd[grepl(clust, siteinfo$station_nm)] <- "well cluster"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="GW")] <- "water table well"
siteinfo$site_tp_cd[which(siteinfo$site_tp_cd=="LK")] <- "surface water (pond)"
# map all sites by type
pal <- colorFactor("Accent", domain = siteinfo$site_tp_cd)

leaflet(siteinfo, options = leafletOptions(zoomControl = FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                 fillOpacity = 1, fillColor = ~pal(site_tp_cd),
                 popup = ~paste(station_nm,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = pal, title = "Site type", opacity = 1, values = ~site_tp_cd, position = "topleft")

\(~\)

2. Specify water quality parameters for data retrieval

# parameters     (nutrients samples are all filtered)
temp <- "00010"  # temperature                 (degC)
sc <- "00095"    # SpC                         (uS/cm)
ph <- "00400"    # pH
do <- "00300"    # DO                          (mg/L)
nh34 <- "00608"  # ammonia and ammonium, as N                (mg/L)
no2 <- "00613"   # nitrite as N                (mg/L)
no3 <- "00618"   # nitrate as N                (mg/L)
no23 <- "00631"  # nitrate + nitrite as N      (mg/L)
tn <- "62854"    # total nitrogen              (mg/L)
po4 <- "00671"   # orthophosphate as P         (mg/L)
delH <- "82082"  # delta H2/H1                (per mil)
delO <- "82085"  # delta O18/O16              (per mil)

pset <- c(temp, sc, ph, do, nh34, no2, no3, no23, tn, po4, delH, delO)

3. Retrieve and summarize data

# RETRIEVAL

# pull water quality data
## option A (to be retired)
data <- readNWISqw(siteNumber = siteinfo$site_no, parameterCd = pset)
b <- data %>% dplyr::select(c(site_no, sample_dt, sample_tm, parm_cd, result_va))

# ## option B (suggested going forward)
# data <- readWQPqw(siteNumber = paste0("USGS-", siteinfo$site_no), parameterCd = pset)
# # select columns and rename
# b <- data %>% dplyr::select(c(MonitoringLocationIdentifier, ActivityStartDate, USGSPCode, CharacteristicName, ResultMeasure.MeasureUnitCode, ResultMeasureValue))
# names(b) <- c("site_no", "sample_dt", "parm_cd", "charname", "units", "result_va")
# b$site_no <- sub("USGS-", "", b$site_no)

# combine site info and qw data
a <- siteinfo 
a$station_nm[grep("SHUBAEL", a$station_nm)] <- "SHUBAEL POND"
alldata <- right_join(a,b) %>%
  rename(date = sample_dt) %>%
  mutate(param = parm_cd, yearmon = format(date, "%Y-%m")) 

# create param names from codes 
alldata$param <- dplyr::recode(alldata$param, `00010` = "temp", `00095` = "SpC", `00400` = "pH", `00300` = "DO", `00608` = "NH34", `00613` = "NO2", `00618` = "NO3", `00631` = "NO23", `62854` = "TN", `00671` = "PO4", `82082` = "delH", `82085` = "delO")

# pull water levels data
gwl <- readNWISgwl(siteNumbers = siteinfo$site_no) 

# join site info
dtw <- gwl %>%
  filter(!is.na(lev_va))
dtw <- right_join(siteinfo %>% dplyr::select(station_nm, site_no, id, dec_lat_va, dec_long_va, alt_va, well_depth_va, weblink), dtw)%>%
  rename(date = lev_dt) %>%
  mutate(yearmon = format(date, "%Y-%m"))
# SUMMARIES

# number of sample depths at each location
nwells <- alldata %>% 
  group_by(id) %>%
  summarise(depths = length(unique(well_depth_va))) 

# number of depths and samples at each location, by parameter
allsum <- alldata %>% 
  group_by(id, param) %>%
  summarise(count = n()) %>%
  pivot_wider(id_cols = id, names_from = param, values_from = count)

allsum <- left_join(allsum, nwells) %>%
  dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
  arrange(depths)

# subset data for a specific sampling event
ym <- "2021-11"
zt <- alldata %>%
  dplyr::filter(yearmon == ym)

# number of depths and samples at each location, by parameter
ztsum <- zt %>% 
  group_by(id, param) %>%
  summarise(count = n()) %>%
  pivot_wider(id_cols = id, names_from = param, values_from = count)

ztsum <- left_join(ztsum, nwells) %>%
  dplyr::select(id, depths, temp, SpC, pH, DO, PO4, NH34:NO3, TN) %>%
  arrange(depths)

# means at each location, single time
zmean <- zt %>%
  group_by(id, parm_cd) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 3))

# means at each location, time series
tsmean <- alldata %>%
  group_by(id, yearmon, param) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=round(mean(result_va, na.rm=TRUE), 3))

# max over sample depths at each location
zmax <- zt %>%
  group_by(id, parm_cd) %>%
  summarise(station_nm=first(id), dec_lat_va=median(dec_lat_va), dec_long_va=median(dec_long_va), weblink=first(weblink), well_depth_va=mean(well_depth_va, na.rm=TRUE), result_va=max(result_va, na.rm=TRUE))

Snapshot from 2021-Nov sampling

# calculate percentages of N in different forms
ftn <- pivot_wider(zt, id_cols = c(id, site_no:dec_long_va, well_depth_va), names_from = param, values_from = result_va)%>%
  filter(!is.na(TN)) %>%
  mutate(percent_organicN = round(100*(1 - ((NH34+NO23)/TN)), 1), percent_NO23 = round(100*NO23/TN, 1))

kable(ftn %>% arrange(id) %>% dplyr::select(station_nm, site_tp_cd, temp, SpC, pH, DO, PO4, NH34, NO3, NO23, TN, percent_organicN, percent_NO23), digits = 2)
station_nm site_tp_cd temp SpC pH DO PO4 NH34 NO3 NO23 TN percent_organicN percent_NO23
MA-A1W 505-0036 well cluster 12.2 143 5.1 9.30 NA 0.02 3.18 3.18 3.23 0.9 98.5
MA-A1W 508-M01-06WT multilevel sampler (MLS) 13.7 105 5.0 3.00 0 0.02 3.28 3.28 3.38 2.4 97.0
MA-A1W 509-0019 water table well 16.8 105 6.3 0.50 0 0.02 0.27 0.28 0.35 14.0 80.3
MA-A1W 514-0028 water table well 11.7 138 4.7 5.50 0 0.02 6.14 6.14 6.30 2.2 97.5
MA-A1W 515-0036 water table well 11.2 214 4.9 7.80 NA 0.02 11.90 11.90 12.40 3.9 96.0
MA-A1W 517-0037 well cluster 11.3 134 4.8 5.50 NA 0.15 3.65 3.65 3.97 4.3 91.9
MA-A1W 520-0048 water table well 10.7 134 5.3 7.80 NA 0.02 6.56 6.56 6.86 4.1 95.6
MA-A1W 521-0041 water table well 10.7 133 5.2 9.60 NA 0.02 1.58 1.58 1.68 4.8 94.0
MA-A1W 522-M01-06WT multilevel sampler (MLS) 16.4 65 5.5 0.23 0 0.02 0.18 0.18 0.23 12.6 78.7
MA-A1W 523-M01-06WT multilevel sampler (MLS) 15.4 260 4.5 0.09 0 1.73 14.10 14.10 16.70 5.2 84.4
MA-A1W 524-0039 well cluster 10.5 151 5.0 3.10 NA 0.02 3.05 3.05 3.24 5.2 94.1
MA-A1W 525-0042 well cluster 11.7 183 4.8 9.10 NA 0.02 8.50 8.50 8.66 1.6 98.2
MA-A1W 526-0034 well cluster 11.9 133 5.8 0.67 NA 0.02 0.42 0.43 0.53 15.8 80.4

\(~\)

\(~\)

4. Map water quality parameters and depth to water

\(~\)

# function to make the maps
mapparam <- function(data, palette, title){
  leaflet(data, options = leafletOptions(zoomControl = TRUE, attributionControl=FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                   fillOpacity = 1, fillColor = ~pal(result_va),
                   popup = ~paste(station_nm,"<br>","Value",result_va,"<br>","<a href='"
                                  , weblink
                                  , "' target='_blank'>"
                                  , "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = palette, title = title, opacity = 1, values = ~result_va, position = "bottomleft")
}

4a. Mean parameter values

z <- zmean

ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1

ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2

ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3

ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4

ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)     
m5 <- mapparam(mapdata, pal, ttl)
#m5

ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6

ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis",  domain = mapdata$result_va)   
m7 <- mapparam(mapdata, pal, ttl)
#m7

ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)    
m8 <- mapparam(mapdata, pal, ttl)
#m8

# map depth to water
pal <- colorNumeric("viridis", domain = dtw$lev_va)
m9 <- leaflet(dtw %>% filter(yearmon == ym), options = leafletOptions(zoomControl = FALSE, attributionControl=FALSE)) %>%
  addProviderTiles('Esri.WorldTopoMap')%>%
  addCircleMarkers(lng = ~dec_long_va, lat = ~dec_lat_va, weight = 1, color = "gray", radius = 6, opacity = 1, 
                 fillOpacity = 1, fillColor = ~pal(lev_va),
                 popup = ~paste(station_nm,"<br>","Value",lev_va,"<br>","<a href='", weblink, "' target='_blank'>", "Link to data</a>")) %>%
  addScaleBar("bottomright") %>%
  addLegend(pal = pal, title = "Depth to water", opacity = 1, values = ~lev_va, position = "bottomleft")

sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)
b = 0

4b. Maximum parameter values

\(~\)

z <- zmax

ttl <- "Temp (°C)"
mapdata <- z %>% filter(parm_cd == temp)
pal <- colorNumeric("viridis", domain = mapdata$result_va)
m1 <- mapparam(mapdata, pal, ttl)
#m1

ttl <- "SpC (uS/cm)"
mapdata <- z %>% filter(parm_cd == sc)
pal <- colorBin("viridis", domain = mapdata$result_va)
m2 <- mapparam(mapdata, pal, ttl)
#m2

ttl = "pH"
mapdata <- z %>% filter(parm_cd == ph)
pal <- colorBin("viridis", domain = mapdata$result_va)
m3 <- mapparam(mapdata, pal, ttl)
#m3

ttl = "DO (mg/L)"
mapdata <- z %>% filter(parm_cd == do)
pal <- colorBin("viridis", domain = mapdata$result_va)
m4 <- mapparam(mapdata, pal, ttl)
#m4

ttl = "Nitrate, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == no3)
pal <- colorBin("viridis", domain = mapdata$result_va)     
m5 <- mapparam(mapdata, pal, ttl)
#m5

ttl = "NH34, as N <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == nh34)
pal <- colorBin("viridis", domain = mapdata$result_va)
m6 <- mapparam(mapdata, pal, ttl)
#m6

ttl = "Total N (mg/L)"
mapdata <- z %>% filter(parm_cd == tn)
pal <- colorBin("viridis",  domain = mapdata$result_va)   
m7 <- mapparam(mapdata, pal, ttl)
#m7

ttl = "Phosphate, as P <br> (mg/L)"
mapdata <- z %>% filter(parm_cd == po4)
pal <- colorBin("viridis", domain = mapdata$result_va)    
m8 <- mapparam(mapdata, pal, ttl)
#m8

sync(m1, m2, m3, m4, m5, m6, m7, m8, m9, ncol = 3, sync.cursor = FALSE)

`\(~\)

b = 0

5. Plot time series

\(~\)

ttl = "Nitrate, as N, (mg/L) in groundwater wells and pond"
p <- ggplot(alldata %>% filter(parm_cd==no3), aes(x = date, y = result_va, group = station_nm)) +
  geom_point(aes(color=id)) +
  geom_line(aes(color=id)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)
ttl = "Mean nitrate, as N, (mg/L) at site"
p <- ggplot(tsmean %>% filter(param=="NO3"), aes(x = yearmon, y = result_va, group = station_nm)) +
  geom_point(aes(color = id)) +
  geom_line(aes(color = id)) +
  theme_bw()+ 
  ggtitle(ttl) 

ggplotly(p)

\(~\)

6. Plot vertical profiles

\(~\)

ttl = "NO3-N (mg/L) vertical profiles"

a <- alldata %>% filter(parm_cd==no3) %>% 
  filter(id %in% c(508, 517, 522, 523, 524, 525)) %>% 
  filter(station_nm != "MA-A1W  517-0036")  # duplicate well @ 517

p <- ggplot(a, aes(x = result_va, y = well_depth_va)) +
  geom_path(aes(color = yearmon))+
  facet_wrap(~id, scales = "free", ncol = 3)+ 
  theme(panel.spacing.x = unit(.25, "cm"), panel.spacing.y = unit(.5, "cm")) +
  #theme_bw() + 
  scale_y_reverse()+ #limits = yrange
  ggtitle(ttl)

ggplotly(p)